function model_output=gfop_soe_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K)

%This function solves for the gains from optimal trade and/or industrial policy,
%starting from the initial equilibrium which is assumed to feature no trade
%or industrial policy.

%CES version with elasticity of substitution rho across sectors

%Function to calculate domestic lambda_hat given wage and employment and
%expenditure shares
function lambda_hat_d=lambda_hat_d_calc(w,lhat,bhat)
    lambda_hat_num =((w*(ones(K,1)./((1+s_k).*lhat.^gamma))).^(-theta)).*lambda;
    lambda_hat_denom=lambda_hat_num+1-lambda;
    lambda_hat_d=lambda_hat_num./lambda_hat_denom.*bhat.*beta;
end

%Function to calculate foreign lambda_hat given wages and employment
function lambda_hat_f=lambda_hat_f_calc(w,lhat)
    lambda_hat_f =(1-t_k).*((w./(((1+s_k).*(1-t_k)).*lhat.^gamma)).^(-theta)).*X;
end

%Function to calculate Revenue given wages and employment and expenditure
%shares
function R=R_calc(w,lhat,bhat)
lambda_hat_d=lambda_hat_d_calc(w,lhat,bhat);
slam=sum(s_k.*lambda_hat_d);
tlamf=sum((t_k.*(1+s_k)-s_k).*((w./(((1+s_k).*(1-t_k)).*lhat.^gamma)).^(-theta)).*X);
R=(-slam*(w+D)+tlamf)/(1+slam);
end

%Function to calculate updated lhat, given wages and initial lhat and
%expenditure shares
function lhat=lhat_calc(w,lhat_in,bhat_in)
    lambda_hat_d=lambda_hat_d_calc(w,lhat_in,bhat_in);
    lambda_hat_f=lambda_hat_f_calc(w,lhat_in);
    R=R_calc(w,lhat_in,bhat_in);
    lhat=((w+D+R)*lambda_hat_d+lambda_hat_f).*((1+s_k)./(w*V_k));   
end

%Function to compute expenditure shares given initial lhat and wages
function bhat=bhat_calc(w,lhat_in)
    lambda_hat_num =((w*(ones(K,1)./((1+s_k).*lhat_in.^gamma))).^(-theta)).*lambda;
    lambda_hat_denom=lambda_hat_num+1-lambda;
    price_hat_term=lambda_hat_denom.^((1-rho)./(-theta));
    bhat=price_hat_term./(sum(price_hat_term.*beta));
end
        

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Algorithm

%Initial values
w=1;
lhat_in=ones(K,1);
bhat_in=ones(K,1);

%Max Iterations
max_lhat_iter=2;   %want to keep this low
max_bhat_iter=2;
max_w_iter=1000000;

%Error tolerance
max_e_w=.000001;
max_e_lhat=.0001;
max_e_bhat=.0001;


%Step sizes
ss_w=.005;
ss_lhat=.005;
ss_bhat=.005;


%%%%%%%%%%%%%
%Algorithm
for i=1:max_w_iter
    for j=1:max_bhat_iter
        for z=1:max_lhat_iter
            lhat= real(lhat_calc(w,lhat_in,bhat_in));
            error_lhat=lhat-lhat_in;
            %zzz=sum(abs(imag(error_lhat))) Diagnostic
            Z_lhat=sum(abs(error_lhat));
            iter_l=z;
            lhat_in=real(lhat_in+ss_lhat*(max(lhat,.01*ones(K,1))-lhat_in));
            lhat=lhat_in;
            if Z_lhat<max_e_lhat
                break;
            end
        end
        bhat=bhat_calc(w,lhat);
        error_bhat=bhat-bhat_in;
        Z_bhat=sum(abs(error_bhat));
            iter_b=j;
            bhat_in=bhat_in+ss_bhat*(bhat-bhat_in);
            bhat=bhat_in;
            if Z_bhat<max_e_bhat
                break;
            end
     end
       %Update wages based on excess demand for labor

       excess_lab_dem=sum(lhat.*V_k)-1;
       Z_w = max(abs(excess_lab_dem));
       w=w+ss_w*(w*excess_lab_dem);
       iterw=i;
       if Z_w<max_e_w
           break;
       end                 
end  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Model Output

%Initializing implicitly calculated variables

R=R_calc(w,lhat,bhat);
lambda_hat_d=lambda_hat_d_calc(w,lhat,bhat);
lambda_hat_f=lambda_hat_f_calc(w,lhat);

%Calculating the change in the price index
lambda_hat_num =((w*(ones(K,1)./((1+s_k).*lhat.^gamma))).^(-theta)).*lambda;
lambda_hat_denom=lambda_hat_num+1-lambda;
price_hat_term=lambda_hat_denom.^((1-rho)./(-theta));

P=(sum(price_hat_term.*beta))^(1/(1-rho));

%Hat change in welfare
welfare=((w+D+R)/(1+D))/P;

%Hat Change in Real income
ri=(w+R)/P;

model_output.lambda_hat_d=lambda_hat_d;
model_output.lambda_hat_f=lambda_hat_f;
model_output.R=R;
model_output.lhat=lhat;
model_output.bhat=bhat;
model_output.w=w;
model_output.P=P;
model_output.welfare=welfare;
model_output.ri=ri;

%Model output for Harberger analysis
model_output.VA_share = lhat.*V_k;   %value added share by sector at the optimum
model_output.lab_real=1-1./lhat;  %percent change in labor allocated to sector k when going from optimum to original equilibrium
model_output.fs_share=(lambda_hat_f)./((1-t_k).*w);   %foreign expenditures on domestic goods, as a share of domestic value added
model_output.fl_reallocation=1-1./(lambda_hat_f./(X.*(1-t_k).*(w./(((1+s_k).*(1-t_k)).*lhat.^gamma))));  %(negative) percent change in effective labor units supplied to foreigners when going from the optimum to original equilibrium

model_output.iterw=iterw;
model_output.Z_w=Z_w
Z_lhat;
Z_bhat;
end